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Abstract 

In this article we extend the results presented in Ref. [Phys. Rev. A 76, 032101 (2007)] to treat 
quantitatively the effects of reservoirs at finite temperature in a bosonic dissipative network: a 
chain of coupled harmonic oscillators whichever its topology, i.e., whichever the way the oscillators 
are coupled together, the strength of their couplings and their natural frequencies. Starting with 
the case where distinct reservoirs are considered, each one coupled to a corresponding oscillator, 
we also analyze the case where a common reservoir is assigned to the whole network. Master 
equations are derived for both situations and both regimes of weak and strong coupling strengths 
between the network oscillators. Solutions of these master equations are presented through the 
normal ordered characteristic function. We also present a technique to estimate the decoherence 
time of network states by computing separately the effects of diffusion and the attenuation of the 
interference terms of the Wigner function. A detailed analysis of the diffusion mechanism is also 
presented through the evolution of the Wigner function. The interesting collective diffusion effects 
are discussed and applied to the analysis of decoherence of a class of network states. Finally, the 
entropy and the entanglement of a pure bipartite system are discussed. 
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I. INTRODUCTION 



The subject of networks of interacting quantum systems has acquired an important role in 
the emerging field of quantum information theory. Since a realistic quantum logic processor 
must ultimately be composed of a large number of interacting quantum systems, it becomes 
mandatory to understand processes as like as perfect state transfer from one to another 
system of the network, and even to compute the fidelity of such a state transfer when the 
action of the environment is taken into account. A significant amount of result has recently 
)een derived on the subject ofperfect state transfer in optical lattices \^ and networks of spin 
2! and harmonic oscillators [sl. Perfect state transfer has also been considered in networks 
of a.b,tra.y topology and coupl.ng configuration fl and even under random fluctuations in 
the couplings of a quantum chains [5|. 

Apart from state transfer, the process of decoherence of a network state has also attracted 
attention and interesting properties of collective damping effects, as the nonadditivity of 
decoherence rates, have been discussed in different contexts as in superconducting qubits 
6], two-atom systems [7|, and chains of dissipative harmonic oscillators ,0,0. still 
regarding the process of collective decoherence, the emergence of decoherence-free subspaces 
(DFSs) has also instigated several interesting results when considering the particular case of 
a composite system interacting with a common reservoir [l2], or the more realistic situation 
where each system interacts with its own reservoir [lOil. We call the attention to the fact that 

nfin 

all Refs. [81, 19|, 110[ envisage such realistic cases of networks where each oscillator interacts 
with its own reservoir, also addressing the particular case where a common reservoir is 
considered. 

To better understand the results in Refs. , [lo|, which are crucial to introduce the 
subject of the present work, we remember that, apart from the distinct reservoirs, the 
network of dissipative harmonic oscillators could present direct an indirect dissipative 
channels. Through the direct channels each oscillator loses excitation to its own reservoir, 
whereas through the indirect channels it loses excitation to all the other reservoirs but to its 
own. When we consider distinct reservoirs for each network oscillator, the indirect dissipative 
channels — intrinsically associated with the nonadditivity of decoherence rates and the 
emergence of DFSs [10] — are significant only in the strong coupling regime where NX^n — 
ue, i.e., the number of network oscillators A^ multiplied by their coupling strengths {Xmn} are 
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about their natural frequencies {ujm}- Therefore, the strong couphng regime, which brings 
together the collective damping effects, depends on the number of network oscillators as 
much as on their coupling strengths. For Markovian white noise reservoirs, however, where 
the spectral densities of the reservoirs are invariant over translation in frequency space, the 
indirect channels becomes null, except for the case = 2 j^. 

In the weak coupling regime where NXmn ^i-, the indirect channels always disappears. 
However, these indirect channels, coming from the strong coupling regime, remains in the 



10|], even assuming a 



case where all network oscillators interacts with a common reservoir 
common Markovian white noise reservoir. This is due to the fact that a common reservoir 
induces an additional correlation between the network oscillators, restoring the indirect 
decay channels. 

Recently, a generalization of Refs. has been presented through a comprehensive 

treatment of networks of dissipative quantum harmonic oscillators, whichever its topology, 
i.e., whichever the way the oscillators are coupled together, the strength of their couplings 
and their natural frequencies HI]. Focusing on the general more realistic scenario where 
each oscillator is coupled to its own reservoir, the case where all the network oscillators 
are coupled to a common reservoir was also addressed. However, after deducing the master 
equation for the case where all the reservoirs are at finite temperatures, all further analysis 
of the dynamics of the network states was restricted to the case where the reservoirs are at 
T = K. Whereas a quantitative analysis of the decoherence and the evolution of the linear 
entropy of representative states of the network were given at K, only a brief qualitative 
analysis of the equilibrium states of the network was presented at finite temperatures. In 
the present manuscript we extend the treatment in Ref. ll| given a detailed analysis of the 
temperature effects on networks of dissipative quantum harmonic oscillators. 

The present extension of Ref. that accounts for the temperature effects coming from 
thermal reservoirs is not only interesting due to its more realistic approach but also from 
the mathematical development here achieved. In fact, we present an alternative approach 
to previous results in the literature [3| regarding the obtainment of the solution of the 
master equation and the estimation of decoherence times through the Wigner distribution 
function. To circumvent noise effects, many of the nowadays experiments demonstrating 
quantum logic operations through atom-field interactions occur in cryogenic environments 
where temperature effects are negligible. In cavity quantum electrodynamics, the setup 
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is cooled to around 0.5 K by a ^He-^He refrigerator to avoid blackbody radiation in the 
High-Q superconducting cavity. Under such a specific condition, the temperature effects on 
the decoherence process are almost negligible. However, when the setup is scaled from one 
single cavity to a network of High-Q cavities, major questions arise due to temperature 
effects. First of all, would the DFSs survive despite the temperature effects? Apart from 
the special class of states composing the DFSs, how the temperature affects other states 
of the network, as for example initial entangled states? Evidently, these questions present 
no obvious answers, even under the assumption that all network cavities are cooled at low 
temperatures. On this regard, we expect that the collective damping effects coming from 
the indirect dissipative channels to play a major role for the answers to the above questions. 

Apart from providing the mathematical treatment of the temperature effects on a network 
of A^ dissipative harmonic oscillators, in the present manuscript we also analyze the role 
played by temperature on the evolution of particular states of the network other than those 
composing DFSs. We reserve the analyses of the emergence of DFSs under temperature 
effects to an specific work [l^ where the mechanism of construction of such privileged states 
is also discussed along with decoherence. 

This paper is organized as following: in Section II we revisit our model of a bosonic 
dissipative network [11] and present the derivation of the master equation governing the 
dynamics of the associated density operator. In Section HI we present the solution of the 
normal ordered characteristic equation obtained from the master equation for the density 
operator of the network. In Section IV we analyze the evolution of two general classes of 
initial states of the network, given by mixtures of coherent and number states, through the 
normal ordered characteristic equation, the Glauber-Sudarshan P-function, and the Wigner 
distribution. A detailed analyses of the diffusion processes is presented in Section V and the 
collective decoherence rates of a family of states of the network is analyzed in Section VI. In 
Section VII we discuss the entropy and the entanglement degree of a pure bipartite system 
and, finally, in Section VIII we present our concluding remarks. 



II. THE MASTER EQUATION OF A BOSONIC DISSIPATIVE NETWORK 



We present here a brief review of the step s for the derivation of the master equation of 
a bosonic network, as developed in Ref. ll|. We start from the general case of a network 
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of interacting oscillators, where each one interacts with each other, from which all other 
topologies can be recovered. As depicted in Fig. 1, we also consider the case where each 
oscillator interacts with its own reservoir due to this more realistic approach for most of 

n 

the physical systems. However, as pointed in Ref. |10|], despite the realistic scenario of the 
case of distinct reservoirs, the case of a common reservoir is more general from the technical 
point of view. In fact, as discussed at the end of this section, the master equation for the 
case of distinct reservoirs can be deduced from the case of a common reservoir. 

We start by considering a general Hamiltonian for a bosonic network, H = Hs + Hr + Hj, 
composed by a set of N coupled oscillators 



N 



m=l 



1 ^ 



(1) 



distinct reservoirs, modeled by a set of A; = 1, . . . , cx) modes. 



N 



Hr = UJmkblnk^mk, (2) 

m=l k 

and the coupling between the network oscillators and their respective reservoirs 



N 

Hi = h'^'^Vmk (&Lam + bmkaC) , 



(3) 



m=l k 



where bl^f^ {hmk) is the creation (annihilation) operator for the kih. bath mode uomk coupled 
to the mth network oscillator uj^ whose creation (annihilation) operator reads a}^ {dm)- The 
coupling strengths between the oscillators are given by the set {Am„}, while those between 
the oscillators and their reservoirs by {Knfc}. We assume, from here on, that i,m,n,i',m', 
and n' run from 1 to A^. 

Before addressing the dissipative process through Hamiltonian ([3]), we focus first on 
Hamiltonian Hs to show how to derive different topologies of a nondissipative network of 
coupled harmonic oscillators. Rewriting Hs in a matrix form 



H, 



h i a\ ■ ■ ■ a]y 



n 
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n 



IN 



\ 



\ 



Hni ■ ■ ■ Ti-NN 
we identify the elements of the matrix Ti = 1-0 a,s 

uOm for m = n 
Xmn for mj^n 



(4) 



(5) 



whose values characterize whichever the network topology, i.e., whichever the way the oscilla- 
tors are coupled together, the set of coupling strengths {Xmn}, and their natural frequencies 

To obtain the master equation of the network we first diagonalize the Hamiltonian Ti, 
(within the physical regime where the normal modes are assumed to be positive) through 
the canonical transformation 

= ^ ] CmnO-ny (6) 
n 

where the coefficients of the mth line of matrix C define the eigenvectors associated to the 
eigenvalues of matrix Ti. With C being an orthogonal matrix, its transposed C" turns 
out to be exactly its inverse C~^, resulting in the commutation relations [Am,^^ = Smn 
and [Am, An] = 0, which enable the Hamiltonian H to be rewritten as a sum H = Hq + V, 
where 



Ho = WmAl^Am + ^ X] ujmkblnkbmk, (7a) 

m m k 

^ = ^Y.Y1 ^^nVmk (hl,A^ + h^uAi) . (7b) 
m,n k 

With the diagonalized Hamiltonian Hq we are ready to introduce the interaction picture, 
defined by the transformation Uo{t) = exp {—iHot/h), where 

V{t) = hJ2 [Omn{t)Ai + OUt)An] , (8) 

m,n 

and Omnit) = C^l^ J2k ^"ik exp [—i {ujmk — '^n) t] bmk- Next, we assume that the interactions 
between the resonators and the reservoirs are weak enough to allow a second-order perturba- 
tion approximation. We also assume a Markovian reservoir such that the density operator of 
the global system can be factorized as psif) ® Pr{0)- Under these assumptions the reduced 
density operator of the network of dissipative coupled resonators satisfy the differential 
equation 

dpsjt) 
dt 

Since for a thermal reservoir {bmkbnk') = (b^mkb\k'^ — 0) have to solve the integrals 
appearing in Eq. ([H]), related to correlation functions of the form 

dr (0.mn{t)0'^„i'ti'^)) = CmnCtm' / dr ''^VmkVm'k' (bmkbm'k'^ 

k,k' 

X exp {-i [{uJrnk - ^n) t - {uJm'k' " ^f) t]} . (lO) 



£ rfr Tr^ [Vit), [V{T),ps{t) ® p«(0)]] . (9) 
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Considering that the reservoir frequencies are very closely spaced to allow a continuum 
summation, we obtain 

dr (Umn{t)U!^,^[T) ) = Ndrnra'C^^Cem 7: 6^ , (11) 

Q \ / Z 

where we have defined the average excitation of the reservoir associated to the mth oscil- 
lator as Um (i^) through the relation {|6j„(z/)6„(z/')) = 27i6mnnm{^)S iy — ^'), apart from the 
damping rates 

7™(^) = fdr r ^ [VMaMf ^e-^(— )(^-*), (12b) 
Jo Jo N-K rimiuj) 

with (7rn{^) being the density of states of the mth reservoir. In the context of Markov 

approximation, where V^(iu„), ami'CUn) and fimi^n) are slowly varying functions around the 

normal modes vUn we can simplify the expressions f lT^ to their usual forms 

7m(t^) = 7m(t^) = [Vra{uj)am{uj)f . (13) 

Back to the Schrodinger picture and to the original field operators am, we finally obtain 
from the steps outlined above, the master equation 
dps{t) i 



m,n 



r + T T 

mn ~ ^ mn f r , \ ^ ran / , \ 
l-'mnPsV') H 7;—'^mnPSV') 



(14) 



where we have defined the damping and the diffusion matrix elements Vmn and T^n in the 
forms 

^mn = ^ X] C'^n7m(cC^^)Cm£, (15a) 
Tmn = ^ X] C'^n7m(^£)^m(tJ7£)Cm£, (15b) 

whereas the Liouville operators accounting for the direct {m = n) and indirect (m 7^ n) 
dissipative channels, are given by 

Cmnpsit) = [anpsit),al^] + [am,ps{t)al] , (16a) 
^mnpsit) = [alpsit),am] + [al^, psit)an] . (16b) 



As mentioned in the Introduction and discussed in Ref. [11|], the oscillators lose excitation 
to their own reservoirs through the direct dissipative channels, whereas through the indirect 



channels they lose excitation to all the other reservoirs but not to their own. Although in 

n 

Ref. [HI] we have obtained the master equation for the general case of reservoirs at finite 
temperatures, all further analysis was carried out for reservoirs at K where the diffusion 
matrix elements T^n are null. Next, considering the case of reservoirs at finite temperatures, 
we must discuss the master equation (1141) under the weak and strong coupling regimes 
between the network oscillators. 



A. Weak coupling regime 

We first remember that the weak coupling regime for a network of coupled oscillator 
is defined by the relation NXmn ^ ^e- (However, if a specific coupling Xmn between two 
oscillators, m and n, fails to satisfy the relation NX^n ^ ^e, the network dynamics is nec- 
essarily described by the strong coupling regime, with some normal-mode frequencies far 
beyond their natural values.) In the weak coupling regime, the interaction between the net- 
work oscillators, described by ^X^m^n -^m" {(4n'^n + (^mdn) /2, could be directly introduced 
into the von Neumann term of the master equation to a good approximation, circumvent- 
ing the necessity to diagonalize the Hamiltonian 7i through a canonical transformation 
= Yln^mnO'n- This is equivalent to approximate the matrix C by a identity matrix I, 
implying that 

Tmn = N'yrn{i^m)5mn, (17a) 
'^mn -^'ym(^^rn)^'m,(^^rn)^rnny (l^b) 

where, evidently, we have also approximated the normal modes by their original natural 
frequencies. Under the above considerations, the master equation f|T^ becomes 



m 



-CmmPs{t) 



H '^mmPsyt) 



(18) 



where, essentially, the indirect dissipative channels disappear, establishing the additivity of 
the decoherence rates. 
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B. Strong coupling regime 



The strong coupling regime means that NX^n ~ ^t, i-e., at least one of the coupling 
{^mn} between two network oscillators must be of the order of any natural frequency {uJm}- 
In this case, the indirect dissipative channels become effective, inducing collective damping 
and diffusion effects that we must investigate. As pondered in the Introduction, would the 
collective effects of nonadditivity of the decay rates and the emergence of DFSs still survive 
despite the temperature effects? 

It must be mentioned that Markovian white noise reservoirs washes out the collective 
damping effects introduced by the strong coupling regime since the spectral densities are 
invariant over translation in frequency space, i.e., jmi^n) = 1m, rendering the same matrix 
elements F^n as in Eq. fll7ap . However, Markovian white noise reservoirs do not washes the 
collective diffusion effects. In fact, only under the additional assumption that nm(ct7„) ^ 
for whatever normal mode vJn-, we recover Eq. fll7bl) for the diffusion matrix T, erasing the 
collective effects completely. Next we discuss the case where the whole network is under the 
action of a common reservoir. 



C. A common reservoir 



When all the network oscillators are coupled to a single reservoir, the master equation. 



derived in Ref. [lij], is similar to that in Eq. (fT^ . replacing the damping and the diffusion 
matrix elements by 

Tmn = N ^ Clnlmn' {'COl) C',], (19a) 

e,n' 

Tmn = N ^ Clnlmn' {'COl) C',], (19b) 

e,n' 

where the damping rates jmni^) and ^ymn (t^) for the case of a single common reservoir are 
given by 

(uj) = dr —V^ (u) Vn {v) a\v)e-'^^-^^\ (20a) 
Jo io N-n 

7„„ (uj) = dr —V^ (u) K (^) a'iu)-^e-''^^--^\ (20b) 
Jo Jo Nn n{uj) 

As mentioned above and discussed in Ref. [lO| , the master equation for the case of distinct 
reservoirs can be deduced from the case of a single common reservoir. The above deduction 



of the master equation (fT^ . where we started from the case of distinct reservoirs, was entirely 
due to its broad apphcation in many physical systems. To demonstrate how to derive the 
case of distinct reservoir from that of a common one, we remember that Vm (i^) gives the 
distribution function of the reservoir modes coupled to the mth oscillator. Therefore, in the 
absence of overlap between the distribution functions, i.e., / c/z/V"m(z/)V„(z/) = for m 7^ n, 
Eqs. ( !20!) reduce to those in Eqs. (fT2!) . In this case, the occurrence of the indirect-decay 
channels follows entirely from the strong coupling between the oscillators, as discussed in the 
subsections presented above. When there is a significant overlap between the distribution 
functions, i.e., / di'Vm{y)Vn{u) 7^ for at least one m ^ n,we get the indirect-decay channels 
even when the network oscillators do not interacts at all. The strength of the damping and 
the diffusion matrix elements being defined by the amount of the overlap, i.e., when the 
overlap between the distributions Vm{i^) and V^(z/) is maximum, the strengths Tmn and T^n 
equals and T^^- 



III. NORMAL ORDERED CHARACTERISTIC FUNCTION 

To analyze the dynamics of the network states for the case where the reservoirs are at 
finite temperatures, we consider the evolution of the (normal ordered) characteristic function, 
derived from the master equation flTll) (suitable for all cases discussed in the previous Section) 

as 



d 



Xi{Vm},t), 



(21) 



where we defined the matrix elements 



T^mn = ^mnf^ + iHr 



(22) 



As noted in Ref. 



the matrix Ti^ is an extension of the free evolution Ti, in Eq. 



which takes into account the dissipative mechanisms of the network. 

Starting with the assumption that Eq. (12T1) admits a solution of the form x{{Vm},t) = 
'^{{Vm})<P{{Vm},t), we obtain two differential equations, one accounting for the dynamic 
process, given by 



d 



-0({r7„},t) = -J] 



d 

Vm [n^^y — + C.C. 



(23) 
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and the other accounting for the stationary solution of the characteristic function, given by 



(24) 



If we perform the substitution — — (7-^^)^ iii^the first differential equation (1231) . it 
turns out to be exactly that appearing in Ref. for the derivation of the solution of 
the Glauber-Sudarshan P-function. Therefore, following the steps outlined in Ref. flll|, the 
solution of Eq. ( l23i) can be written as 



Vr 



.(t) =J2vn (0) Z^:,exp {D^j^y, 



(25) 



£,n 



where we employed the diagonal form of Ti^ following from the transformation D ^ • Ti^ 
JD = fl. Note that writing the solution (l25l) in a matrix form, it becomes 



r^{t) = r/ (0) • D* • exp {-QH) • (D^^)* 
= r/ (0) • exp [- (D • n • D"^)* t] 
= r/(0).exp [- [n'^yt] 



such that 



or, equivalently. 



dr]{t) 
dt 

dr]n{t) _ 
dt 



D\* 



-V it) . (n^) 



(26) 



(27) 



(28) 



representing a system of coupled differential equations which follows from Eq. fl23l) under 



the assumption that (j){{rim},t) = (p{{f]mit)}) , with rjm = rim{0) 

The second differential equation can be solved assuming a general Gaussian form 



^{{Vm}) = exp I -i Er7„n„„< j 

\ m,n / 



(29) 



where the elements of matrix 11 are the coefficients to be determined. Substituting fl29|) into 
Eq. (l2^ and changing conveniently the labels m and n of the involved matrices, we verify 
that the differential equation fl2^ reduces to a matrix equation of the form 



(H^)* • n + n • [n'^y = t + r^, 



(30) 



11 



which is exphcitly written as 

i'^mi) + Umi'H^i = Tmn + '^nm- (31) 



with the superscript T in Eq. (l30l) standing for transposed. It is worth noting that for 
identical reservoirs, where 7m = 7 and so F = F^, we obtain a symmetric dissipative matrix 
H^, i.e., = (H^)^ , making Eq. (150]) the well-known Lyapunov equation. The solution 
of Eq. fl5Ul) . namely the determination of H, can be obtained by converting the matrix 
equation into a system of A^^ algebraic equations, i.e., into a new matrix equation of the 
simplified form A • X = B, with the elements of matrix X being the A^^ unknown variables. 
To this end, it is useful to define the column vector 

vec(n)=(nn Hai ■■■ U^i Hia ■ ■ ■ H^s ■ ■ • Hiat ■ ■ ■ UnnY , (32) 

where the first N elements of vec (11) correspond to the first column of matrix H, whereas 
the next elements correspond to the second column of FI and so on. As so, the matrix 



equation (l30i) can be rewritten into the form 15 1 



[I ® [n^] * + n^(^l]» vec (n) = vec (T + T^) , (33) 

where I is an A x A identity matrix. From the mathematical properties presented in 
Appendix A for the matrix [l ® (H^)* + C?) l] , we verify that the elements of matrix FI 
can be written as 

n«'= E ll^:"" De^D-J^, {DenD-^,)\ (34) 



Q -\-Q* 

m,n,m'n' 



finally leading to the solution of Eq. through Eq. (l2^ . In fact, substituting the expres- 
sion flM|) into the left hand side of Eq. flHTl) . and using the relation D^^ • Ti.^ • D = fl ^ 
•!) = 1) •n, we obtain 

^ ^ \_(P^mi} + '^'mlHn(\ ~ ^ ^ {^m'n' + T^n'm') ^rnn'^nm' — '^mn + '^nmi (35) 

£ m' ,n' 

which is exactly the right hand side of Eq. flHT]) . 

We thus verify that the solution of the characteristic equation is of the form 
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(36) 



Since for t = we get x({^m},0) = ^{{r]m})(p{{'nm},0), such that 0({?7m},O) = 
x{{Vm},0)/(p{{rim}), we end up with the solution of the characteristic function 

xiiVmht) = [xi{Vn.ht = 0)|,,^,^,,,„(,, , (37) 

given in terms of its initial state. An interesting point to be noted is that the dynamics of 
the problem, given by rjmit), takes into account only the dissipative rates Tmn together with 
the free evolution Hamiltonian Ti in Eq. leaving aside the diffusive process associated 
with Tmn- Such a diffusive process, appearing in the ratio {p{{rim}) / ^p{{rimit)}) is however, 
modified, by the dissipative mechanisms. 



IV. DYNAMICS OF THE NETWORK STATES: CHARACTERISTIC FUNC- 
TION, GLAUBER-SUDARSHAN P-FUNCTION, WIGNER DISTRIBUTION AND 
DENSITY OPERATOR 

Starting from two general classes of initial network states, given by mixed superpositions 
of coherent and number states, we next analyze the evolution of such states through the 
characteristic function, the Glauber-Sudarshan P-function, and the Wigner distribution. 
We also compute the network density operator for the case of mixed superposition of Fock 
states. 



A. A mixed superposition of coherent states 

Considering that the initial network state comprehends a mixture of superpositions of 
coherent states like (0))^ = J dvjA (r^) \{Pm i^j)}), the initial density operator becomes 

PsiO) = J2pJ dr,K (r,) j ds,A* (s,) \{(3m (r,)}) {{(3^ is,)}\ , (38) 

where Pj is the probability associated with the state |\1' (0))^. The parameters Vj (sj) rep- 
resent a set of variables defining the probability density function A (r^), and \{l3m i^])}) = 
\(^rn ij j)) stauds for a product of coherent states, where ijj)) represents the state 
associated with the mth network oscillator. In the particular case where A (r^) = 
Ylik^k^ {^3 " '"j'^^)' pure state |\E' (0))^ becomes the discrete superposition 

|vl/(0))^ = 5^A,|{/?^}>^, (39) 

k 
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where we have defined 

(r^^^)]) = Through the definit ion of the time- 

dependent vector elements {ry, t) = J2n ^mnit)Pn i^j) and matrix elements 

Qmn{t) = J2 (-^^^) ^in\ (40a) 

e 

Jmn (t) =Ilran-Yl ^-'n' Q*mm' ) ©nn' (t) , (40b) 
m' ,n' 

we verify, after a rather lengthy calculation, that the evolution of the initial network state 
( l38l) can be described either through the characteristic function 



xiM, t) = J2P^j ds,A is,) J drA* (r,) (r,)} {s,)}) 

X exp J ^ [rimK;^ {rf, t) - ri*^Km {Sf, t)] - VmJmn {t) Vn } , (41) 

I. m m,n J 



either by the Glauber-Sudarshan P-function 



p{{u}, t) = / ^^^^ (^^) / (^^) (^^)> (^^)» 

X exp I -2 5^ (t) - {sy, t)] - ir„ (r,; t)]* I , (42) 

I. m,n ) 

or even by the Wigner distribution function 

w{{U}. t) = I "^'^ (^^) / (^^) (^^)> (^^)» 

det J J J 

X exp I -2 5^ (t) - (s,; t)] [e„ - ir„ (r,; t)]* i . (43) 

I. m,n J 

Note that the difference between the Glauber-Sudarshan P-function and the Wigner 
distribution comes from the time-dependent function associated with the width of the their 
Gaussian function. Consequently, the Wigner function can be obtained from the Glauber- 
Sudarshan P-function through the substitution J — * J = J + I. Whereas the Glauber- 
Sudarshan P-function diverges when there is no diffusion process such that J = (with 
all the reservoirs at K), the width of the Wigner function presents an additional term I 
inhibiting any singularity. 



For the case of K reservoirs 



ll|, the density operator of the network, to be used below, 
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is given by 

B. A mixed superposition of Fock states 

We now assume the initial network state to be a mixture of superposition of Fock states 
1^ — 5^ Cx},...,xj^ \xi, . . . , a; at), where the parameter Xm indicates the number of pho- 

Xl,...,XM 

tons in the mth oscillator while the coefficient Cx},...,xn represents the probability amplitude 
associated with each state . . . ,xn) = \{xm}) composing the whole superposition. The 
initial density operator is thus given by 

PsiO) = E^^E E Kl})*^!!} \i^rn}) {{ym}\ , (45) 

J {Xm}{ym} 

where Pj is the probability associated with the state 1$ (0))^. Since the Fock state \xm), of 
the mth oscillator, can be expanded as a superposition of coherent states of the form 

\xm) = Mm r de^ e-*^-^- \(3^ e^^-> , (46) 

^0 

it is easy to note that the initial state can be obtained by Eq. (1551) . identifying 

A (r,) - A (0,) = E ^{1} n ' (47a) 

{Xm} 

dr^ -> / d9m, (47b) 



Jo 

(3m{rj)^f3me''-, (47c) 

such that we can use the results of the previous subsection to obtain the characteristic func- 
tion, the Glauber-Sudarshan P-function and Wigner distribution for a mixed superposition 
of pure Fock states. Alternatively, such functions may be directly computed from the initial 
state fH5|) . Their expressions are presented in Appendix B, where the density operator for 
a mixed superposition of pure Fock states is also presented. 
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V. TIME EVOLVED DIFFUSION COEFFICIENTS 



To analyze the diffusion mechanism, due to the finite temperature of the reservoirs, we 
start by computing the Wigner distribution associated with the normal-mode oscillators. As 
depicted in Fig. 2, these oscillators, described by Hamiltonian H = Hq + V [Eqs. ( I7a|) and 
(I7bp ]. do not interact with each other, but they do interact with all the reservoirs. We, thus, 
rewrite the Wigner distribution fHSl) in a new coordinate frame {^m}, obtained through the 
diagonalization of matrix J (t). This new framework follows from the rotation 

e = e«U(t); C = v\t)*C, (48) 

where the unitary operation U (t) satisfies U"'^ {t) • J (t) • U (t) = D {t). From this matrix 
relation, we obtain the evolved diffusion coefficients 

'^m (t) = J2 ^In (t) Jnn' it) U^'m (t) (49) 
n,n' 

as the elements of diagonal matrix D (t). In this framework, the rotated Wigner distribution, 
written as 

W{{U}, t) = J2pJ / drjWiiUy, Tj, s„ t), (50) 
are composed by diagonal (r^ = Sj) and off-diagonal (r^ 7^ s^) elements defined by 

W{{Uy. r„ 5„ t) = ^^^A* (r,) A (.,) {{(3^ (r,)} \{(3^ (.,)}) 
det J 

xexpj-^ ^^^^^ ira- Km{sj;t) - ^„ (r^; t) |, (51) 

where K {ry, t)=K {ry, t)»V{t) and K* (r^; t) = {t)»K* {ry, t). The vector K (r^; t) gives 

2 

the excitation intensity of the mth normal-mode oscillator through Km [ry, t) . We stress 
that the larger or smaller values of T>m (t) depend on the network topology (contained within 
the matrix elements Umn), apart from the regime of coupling strengths between the oscillators 
(contained within the matrix elements Jmn {t))- As a particular example of this dependence, 
we consider a degenerate symmetric network, i.e., a degenerate network of oscillators, all 
of them interacting with each other, where {um} = 1^, {^mn} = A, {7™} = 7, {7m} = 7 and 
{nm} = n. In this case, in the strong coupling regime, we obtain the expression 

Jmn (t) = 5mn+'^{l- e-^^*) , (52) 
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and the diffusion coefficients 

P„(t) = | -^-W-^--^^) =1 form^iV, ^^^^ 

[ Jmm it) + {N-l) Jmn (t) = 1 + 2n (l - e"^^*) foi TU = N , 

showing that Jmn (t) 7^ can reduce or enhance the strength of the diffusion coefficients 
(t) associated with the normal-mode oscillators. 



A. Directional and mean diffusion times 

From the above TD diffusion coefficients fl53|) we define the directional diffusion time 

(54) 



1 ^ ^ 

-rPm (t) 



displaying a tendency to a significant spread of the peak — common to all elements (the 
diagonal and off-diagonal) of the Wigner function — associated with the mth normal-mode 
oscillator. Since each normal-mode oscillator defines a direction in the coordinate frame 
{C,m}, we are naturally led to define the mean diffusion time, associated with all the dimen- 
sions of the space, as the average value 



1 1 1 I d ^ , 



(55) 



m diff 

The average diffusion time becomes useful to compute the decoherence time of any network 
state when complemented with the estimated time for a significant decay of the peaks 
associated with the interference terms of the Wigner function (r^ 7^ Sj), to be defined below 

^'S Tint: 

As an illustrative example of the above theory, below we analyze the diffusion coeffi- 
cients T>m (t) for the weak and strong coupling regimes considering the case of a degenerate 
symmetric network. 

1. The weak coupling regime 

In the weak coupling regime, the matrix J{t), already in a diagonal form, is defined 
by the elements Jmn (t) = [1 + 2n (1 — e^'''*)] 6mn, such that U = 1. In this regime, all the 
diffusion coefficients equal to 

Vm (t) =V{t) = l + 2n{l~ e-^') . (56) 
17 



The average diffusion time becomes 

^^^// = (57) 

showing, as expected, that the larger the temperature, the smaller the time required for 
a significant diffusion rate. In this case, the coefficients T>m (t) are mode independent and 
assume a common value, such that the spreads of the peaks associated with the diagonal 
terms of the Wigner function occurs homogeneously in all directions. 



2. The strong coupling regime 

In the strong coupling regime, the elements of matrix J (t) are given by Eq. (152!) and the 
diffusion coefficients by Eq. ( 153|) . showing that only the A^th normal-mode oscillator under- 
goes the diffusion process. For all the normal-mode oscillators but the A^th, the diffusion 
coefficients (t) are counterbalanced by the diffusion rates Jmm (t) and Jmn (t) coming 
from the direct- and indirect-decay channels, respectively. The diffusion coefficients in this 
regime lead to the same mean diffusion time as that in Eq. (1571) . showing that the average 
diffusion effect comes entirely from the temperatures of the reservoirs. As to be demon- 
strated in the next section, this interesting result is not limited to the degenerate symmetric 
topology. 



B. Diffusion and topology 

Starting from Eq. ([55]) and noting that TrS) (t) = Tr J (t) (J (t) = J (t) + I), with the 
elements of matrix J (t) given by Eq. (]40bl) . we obtain the general expression 



-.■^;/ = |TrT, (58) 

applicable to whatever the network topology and the strength coupling regime between the 
oscillators, where 

Tr T = A^ ^ %ni'!^n)nrni'!^n)CnmC~l. (59) 

We note that the information regarding the topology of the network is contained only in the 
product C~l^Cnm which acts as a normalized distribution function {J2m ^nmC^ = 1) when 
computing the average value of the diffusion rate given by Eq. (|58|1 . 
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We identify two general situations where, as in the case of a degenerate symmet- 
ric network, the diffusion mechanism becomes independent of the topology of the net- 
work. The first situation occurs i) when identical reservoirs are assumed, such that 
7m(G7„)ra^(ro„) = 7(tJ7„)n(ti7„) and, consequently, Tr T = N Y,n^i^n)n{^n), making the 
mean diffusion 

T^jf = 2^7(?z7^)n(tz7.^), (60) 

m 

independent of the network topology. The second situation ii) arises from the assumptions 
of Markovian white noise reservoirs and low-temperature regime, where the normal-mode 
frequencies satisfy the relation hzum ^ ksT, ks being the Boltzmann constant. In this 
case we obtain 7m(ro„)nm(tx7„) ^ 7m?^m, such that Tr T = N J^m^rnnm, nm being com- 
puted around the average value of the normal-mode frequencies. The mean diffusion time, 
independent of the network topology, becomes 

m 

Both situation i) and ii) were considered in order to demonstrate that the mean diffusion 
time for both, weak and strong coupling regimes, is the same when considering a degenerate 
symmetric network. For any other situation, apart from i) and ii), the average diffusion 
rate becomes dependent on the network topology, apart from the reservoirs temperatures. 



VI. COLLECTIVE DECOHERENCE RATES 



Since an analysis of decoherence through the density operator of the network is hard 
to derive when temperature effects are present, it becomes appropriate to use the Wigner 
distribution function of the system, instead of the density operator, to estimate the decoher- 
ence time of a family of superposition states which are particular cases of the general state 
given by Eq. (l38l) . This family of states is given by 



\^i,...,Nm=^± 



R 



-a, (3,...,(3) ± 

S N-R-S 



(62) 

where R {S) indicates the number of oscillators in the coherent state a (—a) in the first term 
of the superposition and —a {a) in the second term of the superposition. The remaining 
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N — R — S oscillators are in the coherent state /3. We again stress that we are considering 
a degenerate symmetric network where all the oscillators are indistinguishable. Therefore, 
swapping the states of any two oscillators m and n, we obtain a state which is completely 
equivalent to Eq. flB^ . We also note that when R = 1 and 5 = 0, we obtain from flH^ the 
superposition 

^i,...,iv (0)) = A4 (|«) ± ® \m) , (63) 

where a "Schrodinger cat" -like state is prepared in oscillator 1 while all the remaining oscil- 
lators are prepared in the coherent states (3. 

We start our calculation noting that for a pure two-level state |\&) = a\+) + b\—) , 
whose density matrix is given by p = a*a |+) (+| + b*b |— ) ( — | + a*b |— ) (+| + ab* |+) (— |, 
the ratio of the products between the diagonal and off-diagonal elements equals unity, i.e., 
{a*b) {ab*) / {a*a) (b*b) = 1. For an open system described by a mixed density matrix, 
however, this ratio decrease from unity. Bearing this in mind, we rewrite the Wigner function 
to the discrete case where A (r^) = J2k (^Tj — rj'^'' j , in a form 



(64) 



with its diagonal (r = s) and off-diagonal (r 7^ s) elements given by 



det J 
X JJ^exp 



(65) 



P„ it) 

where r and s (running from 1 to 2) label the product states composing the superposition 
(E2]). 

Now, through the diagonal and off-diagonal elements of the Wigner function, we define 
the ratio 

WrA{L},t)Ws,s{{L},t) 



Wr,s{{U},t)WsA{U},t) 



exp 



(t) 



(66) 



which turns to be independent on the variables {^m} of the Wigner function, as de- 
sired. Moreover, for t = 0, such that 0mn(O) = 5mn and Dm (0) = 1, we obtain 
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Srs(O) = exp (— 1/9^ — . In analogy with the above observation concerning the 
ratio of the products between the diagonal and off-diagonal elements of a pure or mixed 
density matrix, the above defined ratio Srs(t) offers a measure of the decoherence rate 
which follows from the function 



prs (t) = — 



—'rs (t) 



exp 



(67) 



which equals unity for t = 0. The above deduction of the decay function (|6711 can also 
be developed for the general case of an initial continuous superposition state, instead of a 
discrete one. 



A. The equivalence between the decays of the interference terms of both the 
Wigner Function and the density operator: reservoirs at absolute zero 



This subsection is devoted to demonstrate that the measure of the decoherence rate 
offered by Eq. fl67|) is equivalent to the one coming from the interference terms of the density 
operator, which is commonly used for the case of OK reservoirs. In fact, for reservoirs at 
OK, where Dm (t) = 1, it is simple to verify that Eq. (1671) reduces to 



prs (t) = exp 



(68) 



{{Kr^,it)}\{K^it)}) 

where the coefficients {{Pm} \{Pm}) I ({-^m (^)} l{-^m(^)}) those coming from the in- 
terference terms of density operator (l44l) . when considering a discrete case. There- 
fore, considering that decoherence times are usually estimated through the relation 
({/^m} l{/^m}) / ({-^m (''"£')} l{-^m ('^i')}) = i ^he casc of rcscrvoirs at OK, we obtain 



from Eq. (1681) the equivalent relation prs {jd) 
decoherence time through the Wigner function. 



e ^, which gives the estimative of the 
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B. Decay time of the interference terms 



Now we are able to define the time Tint for a significant decay of the peaks associated 
with the interference terms of the Wigner function (r^ 7^ Sj). This is done, by generahzing 
the relation prs (td) = e'^, for the case of reservoirs at finite temperatures, to the equality 



Prs (Tint) = exp 



(69) 



that corresponds to measure the decay of the interference terms of the Wigner function by 
deducting their spreadings, common to all the diagonal and off-diagonal elements, as we can 
see in Eq. fl65|) . In other words, it is similar to analyze the decay of the interference terms in 
a frame where the diagonal terms are frozen. 

C. Decoherence time 

Finally, to define a decoherence time td, which take into account both the diffusion and 
decay of the interference terms, we must consider both the above defined times: the mean 
diffusion time r^iff and the decay time of the interference terms of the Wigner functions 
Tint- We thus define the relation 

(70) 

Td Tdiff Tint 

where T^^jj only becomes relevant for particular initial states whose interference terms of the 
Wigner function are null, as occur, for example in the case = 1, to the coherent state \a), 
or when the excitation of the components of a superposition state is significantly smaller 
than unity. This will become clear in the example to be analyzed below for the degenerate 
symmetric network. In the first case, it is well-known that a coherent state remains as 
such, even under a dissipative process, when considering a reservoir at K. However, when 
considering a reservoir at finite temperature, the decoherence time of a coherent state |a) 
can be estimated through our defined Eqs. fl58|) and fl59l) . 

1. The weak coupling regime 

The Wigner function associated with the pure state (l62l) in the weak coupling regime, is 
obtained from Eq. ([65]) with U = 1 and (t) = V (t) = 1 + 2n (1 - e"^*). Our defined 
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decay function (167|) thus becomes 

Prs {t) = exp [-81?-^ (t) \af {R + S) {1 + 2n) (l - e"^*)] . (71) 

We estimate the decoherence time td of the family of states (!62|) through the equahty 
prs {Tint) = Gxp [— 4P~^ [Tint)]- The obtained resuh for the decay time and so for the deco- 
herence time reads 

_ 1 1 

~ 2|a|S(i?+5)(l + 2n)' ^^^^ 



which recover the results in Ref. lOj for OK reservoirs {n = 0). In Eq. ( 1721) we have disre- 
garded the mean diffusion time r^^^ = 2n7 since we assumed that the excitation {R + S) \af 
is significantly larger than unity. Note that in the case where i? = (5* = 0) or 5* = 
{R = 0), given the initial entangled state V'i,...,Ar (0)^ = Af± (|a, . . . , a, ) ± \—a, . . . , —a)), 
the decoherence time decreases as the number of network oscillators increases. 



For the case of the "Schrodinger cat" -like state in Eq. ( ]63l) . we obtain the result 

'Td ~ Tint = , ,2 , o-N - i'^^) 

Summarizing, the temperature effect decreases the decoherence time when the weak coupling 
regime is considered. 



ipi,...,N (0)), derived from Eq. fl65|) 



2. The strong coupling regime 

From the Wigner function associated with the state 
and using the coefficients (1531) . we obtain in the strong coupling regime 

Prs (t) = exp [-8P^^ (t) (1 + 2n) (l - e"^^*) /N] . (74) 

The estimated decay time Tint of the interference terms of the Wigner functions is established 
through the inequality prs {nnt) = exp {-4N/ [N -1 + Vn {Tint)]} < exp {-41)]^^ [nnt)}, 
such that 

- 2|a|SiV2(l + 2n)' ^^^^ 
showing that the interference terms of the Wigner distribution decay at a fastest rate than in 
the weak coupling regime. For the "Schrodinger cat" -like state, Eq. fl63|) . we obtain exactly 
the result shown in Eq. (!73|) . 
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We finally note that, considering only the usual decay of the interference terms, given 
by prs (td) = exp (—4), the estimation of the decoherence time leads to inconsistent results 
which present negative values apart from singularities. For example, for the "Schrodinger 
cat" -like state in Eq. flB^ . in the particular case = 1, we obtain 

27 [|«|'(l + 2n)-n]' ^^^^ 

which has a singularities at = \a\^ / (l — 2 and becomes negative for \a\^ (1 + 2n) < 
n. Therefore, the procedure adopted in Eq. (!70l) to estimate the decoherence time by 
separating both effects of diffusion and decay of the Wigner function interference terms, 
is in fact more sound than the cruder approach where only the interference effects present 
in the decay function fl67|) are considered. Another example refers to the decoherence of 
a coherent state \a), where the result To ~ 1/2777, computed though the our technique, 
account exactly for the diffusion effect, apart from the decay rate 7, as expected. The usual 
procedure fails to give such an account. 

As mentioned above, the analysis of the emergence of DFSs with the reservoirs at finite 

n 

temperature is addressed in another work p^, where both, collective effects of damping and 
diffusion, are managed together with the network topology to build up desired DFSs. 

VII. COMPUTING THE ENTROPY AND THE ENTANGLEMENT DEGREE 
THROUGH THE WIGNER FUNCTION 

The computation of the density operator of the network for the case of reservoirs at 
finite temperatures becomes a difficult task for the majority of the initial network states. 
Therefore, similarly to our procedure to the analysis of decoherence, we next compute the 
entropy of the network using the Wigner functions as given by 

/oo 
-00 

= l-n^ r Um} W'i{U},t), (77) 
J —00 

where the factor vr^ was introduced to produce a null lower bound for the entropy. Using 
the integral result 

- [ d^r]m exp {amTlla - bmVm - CrnVrnVm) = — eXp ( - ""^"^ ) , (78) 
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and the Wigner function given by Eq. (H3i) . or Eq. (l5Ti) . we obtain the general expression 
S{t) = l- jdrj dr'j ds j ds'A* (r) A* (r') A (s) A (/) {{(3^ (r)} (s)}) 

X ({/?m(r')} |{/3„.(s')})n.s,rv(t), (79) 
which is apphcable to any initial network state, where 

'Prs,r's' (t) 



exp i - ^ [v^ {s, s') v*^ (r, r') 

I. m 



1 



^ Uim {t) ©te (i) Vn {s, s')] Uim {t) Qin {t) Vn (r, r') 



(80) 



i,n / \ i,n 

and t;^ s) = /9m (?") — Pm{s). For the case where dissipation is absent, i.e., 7m(t^) = 
7m(co') = 0, we verify that Vrs,r's' (t) = 1 and, consequently, S = 0. Oppositely, when 
7m(co') 7^ and 7m(^^) 7^ 0, the purity loss follows from the decay of Vrs,r's' (t) which reduces 
to the function p^s {t) ,Eq. fl67|) . that enters in the calculation of the decoherence time, 
under the conditions det J = 1 (T = K), r = r', and s = s'. As expected, the purity loss 
mechanism is intimately related to the decoherence one. 

Focusing on the case when 7m (i^) = lm{^) = 0, the entanglement degree of a bipartite 
system, described by a pure density operator pab ^ A and B standing for two complementary 
sets of network oscillators -, can be computed through the reduced entropy (concurrence) 



C = 1 - TiA [TrB pAsf = 1 - Trs [Tr^ PAsf , 
which is given, through the joint Wigner function W {{C,a} , {(.b} , t), as 



(81) 



C = 1 - TT 



/oo noo 
d'iU} / d'^BW{{U}A^B},t) 

/oo r poo 

-oo J —oo 



(82) 



where A^^ and Nb refer to the numbers of oscillators composing the sets A and B, respec- 
tively. 

When the subsystems A and B are uncorrelated, such that pab = Pa® Pb , the Wigner 
function is factorized as 



and, consequently 
as expected. 



W {{U} , {^b} ,t)=W {{U} , t) W {{^b} , t) , 
C = l-vr-- / d'{^B}Wm^B},t) = 0, 



(83) 
(84) 



25 



VIII. CONCLUDING REMARKS 



In the present work we have analyzed the effects of temperature in a network of dis- 
sipative quantum harmonic oscillators. Starting from a previous work where a general 
treatment of such a bosonic dissipative network was presented ll|] , in the case of reservoirs 
at OK, here we considered reservoirs at finite temperatures. Through the solution obtained 
for the normal-ordered characteristic function, we did compute formal expressions for the 
Glauber-Sudarshan P-function, the Wigner distribution function, and the density operator 
for whichever the initial network state. An important point to be stressed is the relevance 
played by the Wigner function in the present context where the reservoirs are at finite tem- 
perature. In fact, it becomes hard to identify the main features associated with the dynamic 
of the network states through the density operator which results to be an intricate expres- 
sion. Through the Wigner function, however, the diffusion coefficients of the normal-mode 
oscillators are clearly identified as well as the decay of its interference terms. We also showed 
how to compute the entropy and the entanglement degree through the Wigner function. 

We demonstrated that the diffusion coefficients T>m {t) associated with the normal-mode 
oscillators present completely different behaviors in both weak and strong coupling regimes. 
In the former case, where the indirect-decay channels do not take place, the diffusion coeffi- 
cients are entirely related to the dissipative processes of the oscillators to their own reservoirs. 
In this case the collective damping and diffusion effects are dismissible. However, in the later 
case, the diffusion coefficients Vm (t) are counterbalanced by the diffusion rates Jmm (t) and 
Jmn (t) coming from the direct- and indirect-decay channels, respectively. In this case, the 
collective damping and diffusion effects emerges from the fact that all network oscillators 
interact with all the reservoirs due to the strong coupling between each other. In fact, in 
the strong coupling regime, the individual oscillators cannot account for the dynamic of 
the whole network, which must be described through the collective normal-mode oscillators. 
Differently, in the weak coupling regime, the network dynamic follows directly from those 
of the individual oscillators. 

In sum, we have presented an analysis of the mechanisms for handling the diffusion 
coefficients Vm (t) in the strong coupling regime, by manipulating the diffusion rates Jmm (t) 
and Jmn (t) through the nature and the temperature of the reservoirs, apart from the network 
topology. Such approach was explored in Ref. to demonstrate the possibility of the 
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emergence of DFSs in a network of dissipative oscillators even with the reservoirs at finite 
temperatures. 

We have also present a technique to estimate the decoherence time of network states 
which separates effects of diffusion from the decay of the interference terms in the Wigner 
distribution function. Our technique overcomes the difficulties that show up with negative 
values and singularities arising from the usual definition of the decoherence time based 
only on the decay of interference terms. We have computed the decoherence time for some 
particular states of the network, leaving for another work [l^ the analysis of the emergence 
of DFSs under temperature effects. 



APPENDIX A: MATRIX EQUATION 

The solution of an arbitrary matrix equation of the form M • X + X • = P (for an 
unknown X) can be obtained through the solution of the system 

[(I ® M) + (N ® I)] • vec (X) = vec (P) , (Al) 

following from the inverse of [(I ® M) + (N (g) I)], given by 

vec (X) = [(I ® M) + (N O I)]"^ • vec (P) , (A2) 

where the notation vec (P) was defined in Eq.( l32l) . Before computing the elements of the 
inverse matrix (I (g) M) + (N (g) I) = Q, it is useful to observe some important properties of 
Q: 

i) The eigenvalues of matrix Q, defined by obtained through the direct sum of 
the eigenvalues and q of matrices M and N, such that 

eij = ii + ej. (A3) 

ii) The eigenvectors of matrix Q are obtained through the tensor product 

= z>(0 (g, z/O'), (A4) 

where i/^*^ and i>*^*^ describe the eigenvector associated to the eigenvalue and In fact, 
knowing the eigenvalues and eigenvectors of matrices M and N, we can easily verify that 
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'dij defines the desired eigenvector, associated with the eigenvalue Sij, since 



[(I ® M) + (N ® I)] • ( 



) = z>» ® (M«z/(^')) + (N«z>«) 



(A5) 



APPENDIX B: ALTERNATIVE EXPRESSION FOR THE EVOLUTION OF A 
MIXED SUPERPOSITION OF FOCK STATES 



We verify that the evolution of the initial state (H5ll can be characterized, using the same 
definitions (140 al) and fl40bp . through the characteristic function 



X 



iiv^h ^) = E E E^^ (d) * ^{1} (- \ E ^-Jmn it) v: 



n 



J {a;'m}{s/m} 



yi-xp+ji 



Ji 



E 

je=0 



(Bl) 



which leads to the Wigner distribution function given in terms of derivatives as 

[2/nf 



w{{U},t) 



det J 

HE 



rEEEft^ 



.0) \ c^j) 



{y-m} j {Xm} 

1 {Xm}{ym} 



lim 



Qyi-xit+'ik 



X exp 



(B2) 



The above distribution can also be given explicitly in the form 



=1^-^ 5: Eft (c-i'.vc" 



{y,n}l ^{Xm} I 11 



(xe - qeV- 



X 



X 



3 {x,n}{ym} \ t- qf = 

e,e' Ru'=o \ 



J]A,({i?„,„.},{U;t) 



2 J2m.n Jrwn (^) '^■ml{t)^*n 



\ 



R. 



IN- 



exp -2 5^emJ-iC , (B3) 
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where we have defined 



Re,o 

Am{{RiA,{^p};t) 
with 



(^Rn,m—l Rn,m) 



ym—Xm+qm—^l{Rl,m-l^Re,in) 



Rn.m—1 Rn,', 



(B4a) 



, (B4b) 



(B4c) 



Fix) 



(B5) 



1 for X > 
for X < 

As noted in Section IV, we remember that the Glauber-Sudarshan P-function P{{C,m},t) 
can be derived from the Wigner distribution by replacing J^n by J^n. Using such a P- 
function we obtain a compact expression of the evolved density operator associated with the 
initial state (H5l) . given by 



Ps 



X 



3 {Xm}{ym.} 

HE E 



X lim 



e*mQmn{t) j |i„) [ exp | ^ ^ SmJmn (t) E*^ j 

m / J ) \ m,n / 



(B6) 



By defining the parameters 



Rm 


= je + h, 


Sm 


= ii + kg, 


Sen 


= 0, 






6{x) 





if X 7^ 



(B7a) 
(B7b) 
(B7c) 

(B7d) 
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we, alternatively, obtain the explicit form of the density operator 

-EEEp.(C)'^!l,(nf: E 

J {Xm}{ym} 



Ps 



X 



X 



e qi=Oil,jl,kt=0 



{xi - qi)\ 



HE E E 



£' Sii,=0 Kgi,=0 i?^^,=0 



1 



I'-i — on p. 



Sej' — Kfj/V.Kfji] 



X S 



Ri,N — {St 



6 



Hi — xe + qi — ^ {Re'/-i — Re/ 



\ie) {je 
(B8) 

For the case where all the reservoirs are at K, so that Jmn = 0, only the terms with 
Kg// = Se/'^i — Si/r survive in the summation over Ki/i of expression (IBSp . Therefore, at 
K, the density operator flB8l) reduces to the expression 



/ xi 

p^w = EEEft(c^it.)*c!l, HE 



VyJx^- 



E 



3 {Xm}{ym} 

X iJ^iUe} Aki} A)) {^{{Vi- xe + qe} ,t)\ 



nv^^~*Jif— 'n 

1 qe=0 ^ ' ki=0 



(B9) 



as already presented in Ref. 
of product states 



where we have defined, with SgN = 0, the superposition 



i^(te},{M.*)>=0E^^ n E 



xslq,-J2iS, 
\ I' 
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Figure captions 

Fig. 1 Sketch of a dissipative symmetric network of oscillators, where each one interacts 
with each other, apart from its own reservoir. 

Fig. 2 Sketch of a dissipative symmetric network of noninteracting normal-mode 
oscillators, each one interacting with all the reservoirs. 
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Fig. 1 




Fig. 2 



